BRL  TR-3282 


AD-A242  586 


TECHNICAL  REPORT  BRL-TR-3282 


BRL 


A  COMPARATIVE  STUDY  OF  THE 
REACTION  DYNAMICS  OF  A  MODEL  SYSTEM 
USING  DIFFERENT  CRITERIA  IN 
PARAMETERIZING  THE  POTENTIAL  ENERGY  FUNCTION 


BETSY  M.  RICE 
CARY  F.  CHABALOWSKI 
GEORGE  F.  ADAMS 
RICHARD  C  MOWREY 
MICHAEL  J.  PAGE 


OCTOBER  1991 


AI’PROVED  FOR  PUBLIC  RE1J:,\SE;  DISTRIBUTION  IS  UNLIMITED. 


U.S.  ARMY  LABORATORY  COMMAND 


BALLISTIC  RESEARCH  LABORATORY 
ABERDEEN  PROVING  GROUND,  MARYLAND 


91-15881 

li 


NOTICES 


Destroy  this  report  when  it  is  no  longer  needed.  DO  NOT  return  it  to  the  originator. 


Additional  copies  of  this  report  may  be  obtained  from  the  National  Technical  Information  Service, 
U.S.  Department  of  Commerce,  5285  Port  Royal  Road,  Springfield,  VA  22161. 


The  findings  of  this  report  are  not  to  be  construed  as  an  official  Department  of  the  Army  position, 
unless  so  designated  by  other  authorized  documents. 


The  use  of  trade  names  or  manufacturers’  names  in  this  report  does  not  constitute  indorsement 
of  any  commercial  product. 


UNCLASSIFIED 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704  0188 


Public  reponmg  ♦or  tnn  collection  of  information  n  eitimateo  to  average  ’  nouf  oer  resoon^e.  including  tne  time  tor  reviewing  instruaionj,  searcf'ing  e*'Sting  cats  sources, 

gathering  and  maintaining  the  data  needed,  and  como»eting  and  reviewing  the  collection  of  information  Send  comments  regarding  this  burden  estimate  or  any  otr^er  asoect  of  this 
collection  of  information,  including  suggestions  tor  reducing  this  Durden  to  Washington  Headduarters  Services.  Directorate  tor  information  Ooerations  and  Peoorts.  1 2 1 S  Jefferson 
OaviS  Highway.  Suite  t204,  Arlington.  JA  22202-4  302.  and  to  the  Office  of  f^anagement  and  Budget,  ^aoerworn  Reduaion  Projea  (0704-0 '88).  Washington  DC  205(33 


AGENCY  USE  ONLY  (Leave  blank) 


4.  TITLE  AND  SUBTITLE 


2.  REPORT  DATE 

October  1991 


3.  REPORT  TYPE  AND  OATES  COVERED 

Final.  Feb  91  -  Mar  9i 


5.  FUNDING  NUMBERS 


A  (Comparative  Study  of  the  Reaction  Dynamics  of  a  Model  System 
Using  Different  Crrieria  in  Parameterizing  the  Potential  Energy  Function 


6.  AUTHOR(S) 

Betsy  M.  Rice,  Cary  F.  Chabalowski.  George  F.  Adams, 
Richard  C.  Mowrey,  and  Michael  J.  Page 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  AOORESS(ES) 


PR:  1L161102AH43 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


9.  SPONSORING /MONITORING  AGENCY  NAME(S)  AND  AOORESS<ES) 

US  Army  Ballistic  Research  Laboratory 
ATTN:  SLCBR-DD-T 

Aberdeen  Proving  Ground,  MD  21005-5066 


10.  SPONSORING /MONITORING 
AGENCY  REPORT  NUMBER 


BRL-TR-3282 


11.  SUPPLEMENTARY  NOTES 

Richard  C.  Mowrey  and  Michael  J.  Page  are  employees  of  the  Naval  Research  Laboratory,  Washington, 
D.C.  Currently  in  press  for  publication  in  Chemical  Physics  Letters. 


12a.  DISTRIBUTION /AVAILABILITY  STATEMENT 


12b.  DISTRIBUTION  CODE 


Approved  for  public  release:  distribution  is  unlimited. 


13.  ABSTRACT  (Maximum  200  words) 


A  comparative  study  of  the  reaction  dynamics  of  a  model  system  using  a  potential  energy 
function  that  has  been  parameterized  according  to  two  different  criteria  is  presented.  The  molecule 
studied  is  methylene  nitramine  CH^NNOj  and  the  reaction  studied  is  the  N-N  bond  scission  reaction. 
A  mathematical  description  of  the  potential  energy  surface  for  the  system  was  developed  and 
parameters  for  the  individual  terms  in  the  model  were  fitted  to  data  obtained  from  ab  initio 
computations.  The  first  model  was  fitted  to  the  Cartesian  second  derivatives  of  the  energy  for 
equilibrium  methylene  nitramine,  while  the  second  model  was  fitted  using  the  computed  harmonic 
frequencies  for  the  molecule.  Unimolecular  decay  curves  for  the  N-N  bond  scission  reaction  were 
computed  at  six  energies  between  2.9384  and  4.1884  eV.  Not  only  do  the  two  models  disagree  in 
the  computed  reaction  rates,  but  the  rates  predicted  by  Model  1  are  first-order  and  time  independent 
while  those  obtained  with  Model  2  indicate  erratic  and  nonstatist ical  behavior.  This  study  shows  the 
sensitivity  of  dynamics  calculations  to  quantitative  differences  in  potential  energy  surfaces  fitted  to 
different  aspects  of  electronic  structure  calculations. 


14.  SUBJECT  TERMS 


methylene  nitramine;  potential  energy  surface;  molecular  dynamics; 
unirrwlecular  decomposition;  potential  energy;  energy 


17  SECURITY  CLASSIFICATION 
OF  REPORT 

UNCLASSIFIED 


'jSM  /540.0' -/30-SS0Q 


8.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 

UNCLASSIFIED 


SECURITY  CLASSIFICATION 
OF  ABSTRACT 


UNCLASSIFIED 


15.  NUMBER  OF  PAGES 
26 


16.  PRICE  CODE 


20.  LIMITATION  OF  ABSTRACT 


UNCLASSIFIED 


•  t>»o  nv  -‘A 


Intentionally  left  blank. 


II 


TABLE  OF  CONTENTS 


Page 


LIST  OF  FIGURES  .  v 

LIST  OF  TABLES  .  v 

ACKNOWLEDGMENTS .  vii 

1.  INTRODUCTION .  1 

2.  POTENTIAL  ENERGY  SURFACE .  2 

2.1  Form  of  the  Potential  .  2 

2.2  Features  of  the  Two  Models  .  3 

3.  DISCUSSION  .  9 

4.  REFERENCES  .  15 

DISTRIBUTION  LIST .  17 


iii 


■ 


t»r 

■  ^  aJtAAi  ~ 
:  T4ii 


r' 

J  1 1  f  xo  At  1  on. _ _ 


•y _ _ 

Co«M 

jiAb/er 
1,'U. 


;;  I 


Intentionally  left  blank. 


IV 


LIST  OF  FIGURES 


Figure  Page 

1 .  Structure  of  Methylene  Nitramine  at  the  Equilibrium  Geometry  Obtained  at 

the  8-in-8  CASSCF/DZ  Level .  7 

2.  Depictions  of  the  Normal  Modes  of  Vibrations  of  Methylene  Nitramine  for 

Modes  (a)  4,  (b)  11,  and  (c)  12  (See  Table  3) .  8 

3.  Unimolecular  Decay  Curves  for  the  N-N  Bond  Scission  Reaction  of 

Methylene  Nitramine  at  (a)  2.9384  eV  (b)  3.1884  eV,  (c)  3.4384  eV, 

(d)  3.6884  eV,  (e)  3.9384  eV.  and  (f)  4.1 884  eV  Total  Energy .  11 

4.  Plot  of  the  First-Order  Decay  Coefficient  as  a  Function  of  Energy  for 

(a)  Model  1  and  (b)  Model  2  (See  Text) .  12 

LIST  OF  TABLES 

Table  Page 

1 .  Potential  Energy  Functions  and  Parameters .  4 

2.  Structural  Parameters  for  CHjNNOj  at  the  Equilibrium  Geometry  .  5 

3.  Normal  Mode  Frequencies  of  CHjNNOj .  6 

4.  Computed  Rate  Coefficients .  13 


V 


ACKNOWLEDGEMENTS 


The  work  at  the  Naval  Research  Laboratory  (NRL)  was  supported  by  the  Office  of  Naval 
Research  through  NRL;  the  ONR  Physics  Division  (RCM),  Contract  N00014-91-WX-24193; 
and  the  ONR  Mechanics  Division  (MP).  BMR,  CFC,  and  GFA  acknowledge  the  assistance  of 
the  computer  support  staff  of  the  U.S.  Army  Ballistic  Research  Laboratory  (BRL)  and 
equipment  provided  by  the  Productivity  Enhancement  Capital  Investment  Program. 


Intentionally  left  blank. 


1.  INTRODUCTION 


The  desire  to  employ  classical  mechanical  molecular  modeling  and  molecular  dynamics  for 
ever  more  complicated  molecular  systems  has  created  the  need  for  more  general  and 
accurate  descriptions  of  the  potential  energy.  Major  goals  of  both  molecular  mechanics  and 
molecular  dynamics  calculations  are  to  predict  and  understand  the  behavior  of  the  system  at  a 
molecular  level.  The  success  or  failure  in  attaining  these  goals  depends  mainly  on  the 
accuracy  of  the  potential  energy  function  in  describing  the  system.  Often  only  limited  data  are 
available  to  use  in  developing  these  functions,  such  as  vibrational  frequencies  and  activation 
energies  from  experiment,  or  selected  energies  and  geometries  calculated  by  ab  initio 
methods. 

Molecular  mechanics  calculations  employ  force  fields  which  are  usually  tailored  to  describe 
conformational  changes  that  lie  within  a  few  kcal/mol  of  the  equilibrium  structure  and  are 
transferable  among  molecules.  Recently  developed  force  fields,  such  as  MM3  (Allinger,  Yuh, 
and  Lii  1989)  and  CFF89  (Maple  et  al.  1990),  provide  more  general  descriptions  of  materials 
than  did  the  first  generation  of  force  fields.  While  these  force  fields  offer  the  opportunity  to 
model  well  the  properties  of  materials,  there  is  an  increased  requirement  for  molecular 
information  from  which  force  field  parameters  can  be  derived.  In  several  recent  publications, 
force  field  parameters  have  been  derived  by  fitting  to  computed  properties  and  to  the  energy 
derivatives  (first  and  second)  calculated  by  ab  initio  quantum  chemical  techniques  (Maple, 
Dinur,  and  Hagler  1988;  Hagler  et  al.  1989;  Dinur  and  Hagler  1989a,  1989b,  1990).  The  fitting 
techniques  depend  upon  the  application  and  upon  the  desire  to  derive  transferable  force  field 
parameters.  Inclusion  of  the  additional  information  obtained  from  the  ab  initio  calculations  has 
been  reasonably  successful  in  deriving  force  fields  that  better  describe  the  test  molecule  and 
related  molecules. 

In  a  molecular  dynamics  simulation  of  a  molecular  reaction,  by  contrast,  the  energy 
expression  is  generally  designed  for  the  specific  problem,  such  as  a  decomposition  path.  In 
that  case,  the  potential  energy  function  is  designed  to  describe  the  reactant  and  product,  and, 
if  applicable,  the  transition  state  for  the  reaction.  Despite  the  difference  in  focus,  the 
techniques  that  have  provnd  successful  in  the  development  of  force  fields  for  molecular 
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models  should  have  application  in  the  derivation  of  a  molecular  dynamics  potential  energy 
function. 

To  test  this  assertion,  we  construct  potential  energy  expressions  for  the  methylene 
nitramine  molecule,  CHjNNOj,  and  examine  the  effects  of  these  models  on  the  dynamics  of 
unimolecular  decomposition  via  N-N  bond  rupture.  Ab  initio  calculations  for  methylene 
nitramine,  the  transition  state  leading  to  the  elimination  of  MONO  and  the  N-N  bond 
fragmentation  path  have  been  reported  (Mowrey  et.  al  1990).  In  the  present  work,  we 
assumed  a  specific  molecular-mechanics-type  functional  form  for  the  potential  energy  of 
CHjNNOj,  and  the  parameters  in  this  form  were  adjusted  in  two  different  ways.  In  the  first. 
Model  1 ,  the  parameters  were  fit  to  the  matrix  of  second  derivatives  of  the  energy  with  respect 
to  Cartesian  coordinates  computed  in  Mowrey  et  al.  (1990).  The  parameters  for  Model  2  were 
fitted  only  to  the  harmonic  vibrational  frequencies  reported  by  Mowrey  et  al.  (1990).  For  this 
study,  the  properties  were  computed  for  the  equilibrium  structure  as  determined  using  an 
8-in-8  CASSCF  with  a  double  zeta  atomic  basis  set. 

In  Section  2,  we  describe  the  functional  form  of  the  potential  energy  surface  and  comment 
on  the  features.  The  following  section  contains  the  description  of  the  molecular  dynamics 
calculations  for  both  potential  energy  models  and  discusses  the  differences  between  the  two 
sets  of  results. 

2.  POTENTIAL  ENERGY  SURFACE 

2.1  Form  of  the  Potential.  The  function  used  to  describe  methylene  nitramine  is  : 

VT'  I  /Stretch ,  ,  ,  /Stretch  ,  /Stretch  v' 

CH2NNO2  =  I'C/)  +  ^ 

>  ^  E  CT(x,)  ^  [E  CT(x') 

-  f(r^j  , 


-  E  (T,) 

I  /BEND  T->  1  /Torsion  ,  /Torsion  ( 1 ) 

+  ^NN  +  ^  *CNNO  ^  ^HCNN 
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where  all  of  the  interaction  potential  energy  functions  and  parameters  are  given  in  Table  1  and 
F(RfjJ  is  a  function  which  attenuates  interaction  terms  as  N-N  dissociation  occurs.  The  form 
of  F(R^J  is 

1-0.  R.^Re'' 

^(^Nn)  ~  (2) 

where  and  F^'^  are  the  Morse-function  parameters  for  the  N-N  bond  listed  in  Table  1 . 

Some  of  the  terms  in  Equation  1  show  dependence  on  t|  or  T;(the  HCNN  or  CNNO  dihedral 
angles,  respectively).  The  force  constants  and  geometries  differ  for  these  interactions 
depending  on  the  orientation  of  these  dihedral  angles.  The  changes  in  the  two  types  of 
parameters  can  be  described  as  follows; 

p  =  +  (P*^*  -  p“)  [1.0  -  cos(x)]/2  (3) 

where  t  is  the  appropriate  CNNO  or  HCNN  dihedral  angle  and  P  denotes  a  force  constant  or 
geometric  parameter.  The  terms  are  functions  of  the  out-of-plane  wagging  angle,  co. 

The  angle  co  is  given  as 


cos((o)  = 


(R,  X  Rj.) 

^ji  ^jk 


(4) 


and  represents  the  angle  between  a  bond  ji  and  a  vector  perpendicular  to  the  plane  defined 
by  two  bonds  jk  and  jl,  with  co  =  90°  corresponding  to  the  equilibrium  configuration.  The 
potential  energy  term  for  the  wagging  motion  is  the  same  as  for  the  bending  motion. 
Parameters  for  Model  1  were  obtained  by  a  non-linear  least  squares  fit  to  the  Cartesian 
second  energy  derivatives  calculated  by  ab  initio  methods  (Mowrey  et  al.  1990).  Starting  with 
the  parameter  set  for  Model  1 ,  the  parameters  for  Model  2  were  adjusted  by  non-linear  least 
squares  to  fit  the  ab  initio  frequencies  only,  without  regard  to  the  accompanying  eigenvectors. 


2.2  Features  of  the  Two  Models.  The  equilibrium  geometries  and  normal  mode 
frequencies  for  the  two  models  are  shown  for  comparison  with  the  ab  initio  results  in  Tables  2 
and  3,  respectively.  We  have  chosen  the  ab  initio  frequencies  and  eigenvectors  as  our 
standard  for  comparison.  Because  the  parameters  of  the  potential  function  for  Model  2  were 
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Table  1.  Potential  Energy  Functions  and  Parameters 
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Parameter  allowed  to  vary  in  fitting. 
Bond  angle. 

Out-of-plane  wag  angle. 

Same  for  both  models. 


Table  2.  Structural  Parameters  for  CHjNNOg  at  the  Equilibrium  Geometry 


Internal® 

Coordinate 

Ab  initio'’ 

Model  1 

Model  2 

R(CN) 

1.271 

mBH 

1.271 

R(NN) 

1.441 

mSSBm 

1.441 

R(NO)® 

1.303 

BEES 

1.303 

R(NO)'‘ 

1.240 

1.239 

R(CH)® 

1.073 

■Bn 

1.073 

R(CH)' 

1.071 

1.071 

1.071 

0(CNN) 

115.700 

115.700 

115.700 

e(NNO)" 

121.100 

121.100 

121.100 

0(NNO)‘’ 

114.500 

114.500 

114.500 

0(HCN)® 

124.100 

124.100 

124.100 

0(HCN)' 

116.000 

116.000 

116.000 

0(ONO) 

124.400 

124.500 

124.500 

°  Bond  distances  and  angles  are  given  in  units  of  A  and  degrees,  respectively. 
^  8-in-8  CAS3CF/DZ  (Mowrey  et.  al  1990). 

Oxygen  atom  no.  1  (see  Figure  1). 

”  Oxygen  atom  no.  2  (see  Figure  1). 

®  Hydrogen  atom  no.  1  (see  Figure  1). 

'  Hydrogen  atom  no.  2  (see  Figure  1). 


determined  by  fitting  to  the  ab  initio  frequencies,  the  agreement  between  the  two  is  very  close. 
This  is  the  best  fit  to  the  frequencies  that  the  functional  form  would  allow.  The  comparison 
between  the  ab  initio  frequencies  and  those  calculated  from  Model  1  is,  as  expected,  less 
satisfactory.  In  fact,  the  frequencies  calculated  from  Model  1  are  in  significant  disagreement 
with  the  ab  initio  value  for  modes  2,  4,  6,  and  1 1 .  The  eigenvectors  of  Model  1 ,  however, 
fairly  well  resemble  those  of  the  ab  initio  eigenvectors  and  are  easily  assignable  for  all  modes. 
In  contrast,  there  are  6  modes  in  Model  2  (6,  8,  10, 1 1,  12,  and  13)  for  which  the  eigenvectors 
bear  little  resemblance  to  the  corresponding  standard.  Of  these,  modes  8,  10,  11,  and  1 2 
cannot  be  assigned  unambiguously  to  any  of  the  ab  initio  eigenvectors.  Mode  6,  however, 
bears  a  slight  resemblance  to  the  ab  initio  eigenvector  for  mode  5.  To  illustrate,  the 
eigenvectors  corresponding  to  modes  4,  1 1 ,  and  12  calculated  by  the  2  models  are  shown  for 
comparison  with  the  standards  in  Figure  2.  Figure  2a  shows  the  worst  agreement  in 
frequency  between  Model  1  and  the  standard.  Inspection  of  this  figure  shows  that  the 
eigenvectors  of  both  of  the  models  are  in  relatively  good  agreement  with  the  standard. 
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Table  3.  Normal  Mode  Frequencies  of  CHjNNO® 


Mode  No. 

Ab  initio*’ 

Model  1' 

Model  2'’ 

1 

100 

101 

101 

2 

391 

274 

391 

3 

561 

560 

560 

4 

621 

441 

619 

5 

665 

610 

666 

6 

872 

732 

868 

7 

896 

894 

894 

8 

1220 

1206 

1220 

9 

1261 

1263 

1263 

10 

1311 

1312 

1311 

11 

1580 

1447 

1580 

12 

1612 

1637 

1613 

13 

1839 

1839 

1840 

14 

3361 

3362 

3356 

15 

3483 

3492 

3487 

ZPE®  (eV) 

1 .2259 

1.1884 

1.2255 

®  Frequencies  given  in  cm  '. 

*’  8-in-8  CASSCF/DZ  (Mowrey  et  al.  1990). 

'  Eigenvectors  corresponding  to  these  frequencies  matched  to  ab  initio  eigenvectors. 
“  Frequencies  matched  to  ab  initio  frequencies  regardless  of  eigenvector  form. 

®  Zero  point  vibrational  energy. 


Figure  2b  again  demonstrates  a  case  where  the  frequency  of  the  mode  predicted  by 
Model  1  disagrees  with  the  standard,  but  the  eigenvector  is  clearly  assignable.  In  this  case, 
however,  the  eigenvector  of  Model  2  does  not  resemble  the  standard  in  spite  of  its  excellent 
agreement  in  frequency.  In  the  last  case,  illustrated  in  Figure  2c,  the  frequencies  from  both 
models  compare  well  with  the  standard,  but  the  eigenvector  of  Model  2  is  unassignable  to  this 
mode.  The  differences  in  vibrational  frequencies  between  Model  1  and  the  standard  merely 
reflect  the  inability  to  perfectly  fit  the  matrix  of  second-derivatives  with  respect  to  energy  for 
Model  1  which  indicates  deficiencies  in  the  form  of  the  potential  energy  function.  On  the  other 
hand,  the  inability  to  reproduce  the  eigenvectors  from  Model  2  is  due  not  only  to  deficiencies 
in  the  functional  form,  but  also  to  the  fact  that  the  parameters  for  Model  2  were  fit  to  less 
potential  energy  surface  information,  specifically  concerning  the  shape  of  the  surface. 
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Ab  initio  results  are  shown  in  the  first  column  of  the  figures,  while  the  results  for  Models  1  and  2  are  shown  in  the 
second  and  third  columns,  respectively.  The  harmonic  frequencies  of  the  modes  are  listed  in  parentheses. 


Figure  2.  Depictions  of  the  Normal  Modes  of  Vibrations  of  Methylene  Nitramine  for  Modes 
(a)  4,  (b)  11,  and  (c)  12  (See  Table  3). 
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In  almost  all  instances,  there  will  be  deficiencies  in  the  form  of  the  potential  energy 
function;  therefore,  it  will  be  difficult  to  attain  good  fits.  In  addition  to  defining  a  suitably 
flexible  function,  a  modeler  is  also  faced  with  the  issue  of  whether  it  is  better  to  correctly 
describe  the  eigenvectors  at  the  sacrifice  of  agreement  in  frequencies  or  to  fit  to  the 
frequencies  only,  thereby  ignoring  the  detailed  motion  of  the  atoms.  For  this  system,  the 
assumption  that  Model  2  accurately  describes  the  system  under  study  based  only  on 
agreement  with  the  harmonic  frequencies  is  deceptive,  considering  how  poorly  some  of  the 
eigenvectors  compare  with  the  ab  initio  values.  Although  differences  in  the  dynamics  occur 
near  the  equilibrium  geometry  due  to  the  differences  in  the  fitting  criteria  for  the  model 
(i.e.,  the  normal  modes),  we  wish  to  further  explore  how  the  reaction  dynamics  differ 
depending  on  the  different  criteria  used  in  fitting  the  functions.  To  investigate  this,  we 
performed  a  series  of  molecular  dynamics  calculations  on  the  two  models;  specifically,  we 
investigated  the  unimolecular  decomposition  of  CHgNNOj. 

3.  DISCUSSION 

Although  there  are  experimental  evidence  (Zhao,  Hintsa,  and  Lee  1988)  and  theoretical 
predictions  (Mowrey  et.  al  1990)  that  unimolecular  decay  of  CHjNNOj  is  more  likely  to  occur 
by  concerted  dissociation  pathways  than  through  N-N  bond  scission,  for  our  purposes  our 
models  only  describe  bond  scission  reactions.  The  N-N  bond  is  the  weakest  bond  in  this 
system  by  2.4  eV,  and  this  was  the  only  reaction  we  observed. 

Ensembles  of  trajectories  at  six  energies  were  integrated  using  a  variable  step  size 
Adams-Moulton  fourth-order  predictor  corrector  integrator  (Miller  and  George  1972).  Relative 
error  tolerance  was  set  at  10  ^  Before  each  ensemble  of  trajectories  is  integrated,  the  system 
(including  the  zero  point  energy)  is  equipartitioned  among  the  normal  modes  in  the  form  of 
kinetic  energy.  A  warm-up  trajectory  of  approximately  0.03  ps  is  performed,  and  a  Markov 
walk  of  75,000  steps  is  taken  to  randomize  the  energy  of  the  system.  A  trajectory  is  then 
integrated  until  the  N-N  bond  exceeds  6.0  A  or  until  the  trajectory  integration  exceeds  a 
maximum  time  of  30  ps.  A  sequence  of  1,500  Markov  moves  is  taken  from  the  starting  point 
of  the  previous  trajectory  and  the  integration/Markov  walk  pattern  is  repeated  until  1 ,000 
trajectories  have  been  integrated.  According  to  the  standard  Metropolis  Monte  Carlo  sampling 
procedure  (Metropolis  et.  al  1953;  Raff  and  Thompson  1985;  Brady,  Doll,  and 
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Thompson  1981),  if  a  trial  move  taken  from  the  starting  point  of  a  previous  trajectory  j  is 
rejected,  then  the  result  of  the  previous  trajectory  is  counted  again  and  recorded  as  the  j  +  1 
contribution  to  the  ensemble  average.  The  present  Metropolis  procedure  was  designed  such 
that  approximately  50%  of  all  attempted  moves  were  rejected.  Because  of  this  rejection  ratio, 
results  from  approximately  50%  of  the  1 ,000  trajectories  integrated  were  recounted  and 
included  in  the  decay  curves.  In  the  Markov  walk,  the  system  is  not  allowed  to  exceed  an 
N-N  bond  distance  of  3.75  A.  Each  trajectory  included  in  the  decay  curves  was  required  to 
undergo  at  least  one  N-N  bond  vibration  before  dissociation;  if  it  did  not,  it  was  not  included. 
The  time  of  decomposition  was  marked  at  the  last  inner  turning  point  of  the  N-N  bond 
vibration. 


Figure  3  shows  unimolecular  decay  curves  of  ln(P)  vs.  time  where  P  is  the  fraction  of 
undissociated  trajectories  at  time  f  for  Models  1  and  2  at  2.9384,  3.1884,  3.4384,  3.6884, 
3.9384,  and  4.1884  eV.  First-order  rate  coefficients  were  extracted  from  the  straight  line 
portions  of  the  curves  and  are  listed  in  Table  4.  In  all  of  the  figures,  the  decay  curves  of 
Model  1  are  well  fit  to  a  straight  line,  indicating  the  rate  in  this  system  is  first-order  and  time 
independent.  Also,  as  expected,  the  rate  coefficients  show  a  monotonic  increase  with 
increasing  energy.  Another  indication  of  the  well-behaved  nature  of  this  model  is  that  the 
rates  can  be  described  by  the  RRK  statistical  model  (Robinson  and  Holbrook  1972) 


k{E)  =  A 


E  -  E, 


(4) 


where  A  =  5.1  ps  \  E^  =  2.7  eV  and  s  =  1.3.  The  calculated  rates  from  Model  1  and  those 
predicted  from  Equation  4  are  shown  in  Figure  4a. 


The  decay  curves  for  Model  2  are  not  well  behaved;  in  Figures  2b,  2d,  and  2e  there  are  at 
least  two  regions  that  are  linear  but  differ  in  slope.  The  slopes  for  the  initial  regions  of  these 
three  curves  are  significantly  sleeper  than  for  the  region  of  longer  lifetimes,  indicating  that 
there  is  a  faster  reaction  rate  for  the  early  decomposition  events  and  a  slower  rate  for  longer 
reaction  times.  Like  Model  1 ,  the  three  other  decay  curves  (Figures  2a,  2c,  and  2f)  can  be 
reasonably  fit  to  a  straight  line.  We  attempted  to  fit  the  data  from  Model  2  to  Equation  4  in 
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Figure  3.  Unimolecular  Decay  Curves  for  the  N-N  Bond  Scission  Reaction  of  Methylene  Nitramine  at  (a)  2.9384,  (b)  3.1884  eV 
(c)  3.4384  eV.  id)  3.6884.  (e)  3.9384  eV.  and  (f)  4.1884  eV  Total  Energy. 


3.0  3.3  3.6  3.9 

ENERGY  (eV) 


The  filled  circles  are  the  rates  extracted  from  the 
these  results  to  the  RRK  expression,  Equation  4. 


trajectory  results  and  the 


curves  are  the  results  of  the  best  fit 


of 


Figure  4. 


P.lpt  of  the  First-Order  Decay  Coefficient  as 
Lb)  Model  2  (See  Text)  ' 


a  Function  Energy  for  fa)  Model  1  anH 
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Table  4.  Computed  Rate  Coefficients' 


Energy  (eV)*’ 

k  (ps 

Model  1 

Model  2" 

Fast 

Slow 

2.94 

1.64 

3.96 

3.19 

2.54 

1.46 

0.142 

3.44 

2.98 

4.36 

3.69 

3.23 

3.71 

1.310 

3.94 

3.48 

4.80 

0.729 

4.19 

3.57 

5.38 

The  rate  coefficients  k  are  extracted  from  the  least  squares  fits  of  the  decay  curves  to  the  model  of  first  order 
reaction,  ln(p)  -  -kt  (Figure  3). 

Includes  the  zero  pxDint  energy. 

See  text. 


the  same  manner  as  for  Model  1  using  the  rate  coefficients  extracted  from  the  curves  in 
Figures  2a,  2c,  and  2f  and  the  "fast"  region  of  the  curves  in  Figures  2b,  2d,  and  2e.  Because 
Figures  2a,  2c,  and  2f  did  not  unambiguously  exhibit  a  "slow"  region  of  the  curves,  we  did  not 
attempt  to  fit  the  "slow"  rates  to  Equation  4.  Comparison  of  the  best  fit  of  Equation  4  to  the 
rates  extracted  from  Model  2  are  shown  in  Figure  4b;  the  best  fit  parameters  for  this  equation 
are  11.9  ps  \  1.178  eV  and  2.97  for  A,  E^,  and  s,  respectively.  As  shown  strikingly  in 
Figure  4b  and  Table  4,  the  rates  from  Model  2  have  an  erratic  energy  dependence,  and  are 
poorly  described  by  the  RRK  theory.  These,  coupled  with  bimodal  behavior  of  the  decay 
curves  in  Figures  2b,  2d,  and  2e,  are  indicative  of  nonstatistical  behavior,  possibly  due  to 
bottlenecks  in  the  system  (Hase  1 976). 

Because  CH^NNOj  has  never  been  isolated  and  unimolecular  decay  experiments  have 
never  been  done  on  this  molecule,  particularly  this  channel,  we  cannot  state  unequivocally 
which  model  more  accurately  represents  the  true  system.  We  have  shown,  however,  that  with 
little  effort,  nonstatistical  behavior  can  be  incorporated  into  a  model,  and  at  the  very  least,  we 
have  shown  that  the  rates  of  decomposition  differ  between  the  two  models  at  all  energies 
studied.  We  wish  to  stress  that  Model  2,  which  predicts  nonstatistical  behavior,  is  fit  to  the 
observable,  the  normal  mode  frequencies.  If  this  nonstatistical  behavior  is  only  an  artifact  of 
the  fitting  procedure,  one  could  easily  be  led  to  erroneous  conclusions  about  a  given  system 
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which  exhibits  very  interesting  dynamical  behavior.  The  erratic  nature  of  the  dynamics  results 
obtained  from  Model  2  leads  one  to  question  the  usefulness  for  prediction  and  understanding 
the  system  as  described  by  Model  2. 
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